clc;
% Read Excel LP data
str='q.xlsx';
% str='t.xlsx';
% str='AFIRO.xlsx';
c=xlsread(str, 'c');
ALEqual=xlsread(str, 'A(<=)');
bLEqual=xlsread(str, 'b(<=)');
AGEqual=xlsread(str, 'A(>=)');
bGEqual=xlsread(str, 'b(>=)');
Aeq=xlsread(str, 'Aeq(=)');
beq=xlsread(str, 'beq(=)');
lb=xlsread(str, 'lb');

tic;
% Get no. of constraints and variables
global cRow;
[cRow cCol]= size(c);                
[ALEqualRow ALEqualCol]= size(ALEqual);                 % <=
[bLEqualRow bLEqualCol]= size(bLEqual);                 % <=
[AGEqualRow AGEqualCol]= size(AGEqual);                 % >=
[AeqRow AeqCol]= size(Aeq);                             % =

constraint = ALEqualRow + AGEqualRow + AeqRow;        % m
variable = max(ALEqualCol,max(AGEqualCol,AeqCol));    % n

noIneqConstraint = ALEqualRow + AGEqualRow;

slackVar = noIneqConstraint + variable;        % slack variable

% Transform import LP to standard form Ax <= b, Bx = 0
% Add zeros for c, B, Aeq, b
global nC nAeq nbeq;
nC = c;
nALEqual = ALEqual;
nAGEqual = AGEqual;
nAeq = Aeq;
nbeq = beq;

% New beq - combine sum of inequality and equality constraints
nbeq = [beq;bLEqual];
nbeq = [nbeq;bGEqual];

% Find no. of variables for displaying result
global cVariable;
cVariable = [];
n = 1;
for ii=1:1:cRow
    if c(ii) ~= 0
        cVariable(n,1) = ii;
        n = n + 1;
    end    
end

% New c
nC(slackVar,1) = 0;

% New Aeq
% (=)constraint
if (isempty(nALEqual)==0) || (isempty(nAGEqual)==0)
    if isempty(nAeq) == 0
        nAeq(AeqRow,slackVar)= 0;
    end
end

% (<=)constraint
if isempty(nALEqual) == 0
    nALEqual(ALEqualRow,slackVar)= 0;
    % add inequality constraints to new A
    nAeq=[nAeq;nALEqual];
end

% (>=)constraint
if isempty(nAGEqual) == 0
    nAGEqual(AGEqualRow,slackVar)= 0;
    % add inequality constraints to new A
    nAeq=[nAeq;nAGEqual];
end

% index of A (constraints) row
if isempty(nAeq)
    rowOfAeq = 1;
else
    rowOfAeq = 1 + AeqRow;
end

% index of A (constraints) column (slack variable)
colOfAeq = variable + 1;

% Aeq - add slack variable for inequality constraints forchanging to equality constraint
% New A & b - add inequality constraints for those constraints

global nA nb;
nA=[];
nb=[];
[lbRow lbCol] = size(lb);
nARow = 1;

% (=)constraint
% nA = zeros(colOfAeq-1,1);

% New A (lower bound from excel)
for ii=1:1:lbRow
    nA(nARow,nARow) = -1;
    nARow = nARow + 1;
end

% New b (lower bound from excel)
% nb = zeros(colOfAeq-1,1);
nb = lb;


% (<=)constraint
if isempty(nALEqual) == 0
    for ii=rowOfAeq:1:AeqRow + ALEqualRow
        nAeq(rowOfAeq,colOfAeq) = 1;
        nA(nARow,nARow) = -1;
        nb=[nb;0];
        rowOfAeq = rowOfAeq + 1;
        colOfAeq = colOfAeq + 1;
        nARow = nARow + 1;
    end
end

%(>=)constraint
if isempty(nAGEqual) == 0
    for ii=rowOfAeq:1:AeqRow + ALEqualRow + AGEqualRow
        nAeq(rowOfAeq,colOfAeq) = -1;
        nA(nARow,nARow) = -1;
        nb=[nb;0];       
        rowOfAeq = rowOfAeq + 1;
        colOfAeq = colOfAeq + 1;
        nARow = nARow + 1;
    end
end

toc

tic;
% [x,fval,exitflag,output,lambda]=linprog(nC,nA,nb,nAeq,nbeq)
[x2,fval,exitflag,output,lambda]=linprog(nC,nA,nb,nAeq,nbeq)

toc

% xlswrite('ExportA.xlsx', nC, 'nC');
% xlswrite('ExportA.xlsx', nA, 'nA');
% xlswrite('ExportA.xlsx', nb, 'nb)');
% xlswrite('ExportA.xlsx', nAeq, 'nAeq');
% xlswrite('ExportA.xlsx', nbeq, 'nbeq');

n=1;
nX2=[];
temp=0;

for ii=1:1:cRow
    if n <= length(cVariable)
        if ii == cVariable(n)
%             if x2(ii) < 0.0001
%                 temp = 0;
%             else
%                 temp = (x2(ii)*100)/100;
%             end
            temp = x2(ii);
            nX2{n,1} = ['x' num2str(ii) ' = ' num2str(temp)];
            n = n + 1;
        end
    end
end

nX2

% xlswrite('Da.xlsx', nX2, 'x');
% ======================
tic;
%------------------------------M, r, y
[nrowAeq, ncolAeq] = size(nAeq);

% check Aeq has full row rank
if rank(nAeq)==nrowAeq
    fullRowRank=true;
end

if fullRowRank
    % generate a random n x n non-singular matrix M
    M = randi([0,5],ncolAeq);
    
    while det(M)==0
        M = randi([0,5],ncolAeq);
    end
    
    % generate a random n x 1 vector r
    r = randi([0,5],ncolAeq,1);
    
    while nbeq + nAeq * r == 0
        r = randi([0,5],ncolAeq,1);
    end
    
    % generate a random integer between 1 and 50
    gamma = round(rand(1)*49)+1;
    
    cP = gamma * M' * nC;           % c' = gM'c
    
    % -----------------------Q, A', b'
    % generate a m x m non-singular matrix Q (deteminant =/= 0)
    Q = randi([0,5],nrowAeq);
    
    while det(Q)==0
        Q = randi([0,5],nrowAeq);
    end
    
    AeqP = Q * nAeq * M;                    % A' = QAM
    beqP = Q * (nbeq + (nAeq * r));         % b' = Q(b + Ar)
        
    %-----------------------------B'
    [nrowAeqP, ncolAeqP] = size(AeqP);
    nonZero = nrowAeqP;         % find the non-zero of nbeq from end of the vector
    allZero = true;             % determine nbeq is all zero
    
    % find the non-zero of nbeq from end of the vector
    while nonZero ~= 0
        if beqP(nonZero) ~= 0
            allZero = false;            % found & exit loop
            break;
        end
        nonZero = nonZero - 1;
    end
    
    % generate a random n x m lambda
    Lambda =randi([0,5],ncolAeqP,nrowAeqP);
    sum = 0;      % sum of each row
    
    Br = nA * r;              % Br
    % change content of Lambda in order to cause Lambda * b' = Br
    if allZero == false
        if nonZero == 1
            % change first column of Lambda to zero
            for ii=1:1:ncolAeqP                  % row
                if nonZero == 1
                    Lambda(ii, 1) = 0;
                end
            end
        else
            for ii=1:1:ncolAeqP                  % row
                for jj=1:1:(nonZero-1)           % col
                    sum = sum + (Lambda(ii,jj) * beqP(jj));
                end
                Lambda(ii, nonZero) = (Br(ii) - sum) / beqP(nonZero);
                sum = 0;
            end
        end
    end    
 
    AP = (nA * M) - (Lambda * AeqP);         % B' = BM - (Lambda)QAM => BM-LA' =>(A' = AeqP)
    bP = (nA * r) - (Lambda * beqP);         % 0 = Br - Lambda(Q(b+Ar)) =>(Q(b+Ar) = beqP)

    % display transform Time
    toc
    
    options = mskoptimset('MaxIter',4000000, 'Write', 'prob.opf');
    
    tic;
    [y,fval,exitflag] = linprog(cP,AP,bP,AeqP,beqP,[],[],[],options);
    toc
        
    % x2 = round (((M * y) - r) * 100 ) / 100;
    x = (M * y) - r;
    
    n=1;
    nX=[];
    temp=0;
    
    for ii=1:1:cRow
        if n <= length(cVariable)
            if ii == cVariable(n)
%                 if x(ii) < 0.0001
%                     temp = 0;
%                 else
%                 temp = (x(ii)*100)/100;
                temp=x(ii);
%                 end
                nX{n,1} = ['x' num2str(ii) ' = ' num2str(temp)];
                n = n + 1;
            end
        end
    end
    nX
end